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ABSTRACT 

We compare the substructure evolution in pure dark matter (DM) halos with those in the presence 
of baryons, hereafter PDM and BDM models. The prime halos have been analyzed in the previous 
work, Romano-Diaz et al. Models have been evolved from identical initial conditions which have been 
constructed by means of the Constrained Realization method. The BDM model includes star formation 
and feedback from stellar evolution onto the gas. A comprehensive catalog of subhalo populations has 
been compiled and individual and statistical properties of subhalos analyzed, including their orbital 
differences. We find that subhalo population mass functions in PDM and BDM are consistent with 
a single power law, M^ h , for each of the models in the mass range of ~ 2 x 10 s M — 2 x 10 11 M Q . 
However, we detect a nonnegligible shift between these functions, the time-averaged a ~ —0.86 for the 
PDM and —0.98 for the BDM models. Overall, a appears to be a nearly constant with variations of 
±15%. Second, we find that the radial mass distribution of subhalo populations can be approximated 
by a power law, i? 7abh with a steepening that occurs at the radius of a maximal circular velocity, i? V max, 
in the prime halos. Here we find that the 7 s bh ^ — 1.5 for the PDM and -1 for the BDM models, when 
averaged over time inside R vmax . The slope is steeper outside this region and approaches —3. We 
detect little spatial bias (less than 10%) between the subhalo populations and the DM distribution of 
the main halos. Also, the subhalo population exhibits much less triaxiality in the presence of baryons, 
in tandem with the shape of the prime halo. Finally, we find that, counter-intuitively, the BDM 
population is depleted at a faster rate than the PDM one within the central 30 kpc of the prime 
halo. The reason for this is that although the baryons provide a substantial glue to the subhalos, 
the main halo exhibits the same trend. This assures a more efficient tidal disruption of the BDM 
subhalo population. However, this effect can be reversed for a more efficient feedback from stellar 
evolution and the central supermassive black holes, which will expel baryons from the center and 
decrease the central concentration of the prime halo. We compare our results with via Lactea and 
Aquarius simulations and other published results. 

Subject headings: cosmology: dark matter — galaxies: evolution — galaxies: formation — galaxies: 
halos — galaxies: interactions — galaxies: kinematics and dynamics 



1. INTRODUCTION 

The high-redshift Universe is characterized by a uni- 
form mixture of dark matter (DM) and baryons (e.g., 
Spergel et al. 2007). The process of galaxy formation 
is expected to lead essentially to a (partial) separation 
between the DM and baryons, due to dissipative pro- 
cesses in the latter. As a result, the central regions of 
many galaxies should be dominated by the baryons and 
by their dynamics — this is supported by observations 
of various aspects of galactic dynamics (e.g., Flores & 
Primack 1994; de Blok & Bosma 2002; Sand et al. 2004; 
Gentile, Tonini & Salucci 2007; de Blok et al. 2008). 
On the other hand and as a direct consequence of the 
process of a hierarchical buildup of structure in the Uni- 
verse, a large amount of substructure, i.e., subhalos, pen- 
etrating the prime halos is expected to be present (e.g., 
Klypin et al. 1999; Moore et al. 1999). Moreover, the 
baryons are expected to modify the DM structure on sub- 
galactic scales, both in the prime halos and the subhalos, 
although the extent of this change is far from being clear. 
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Such differences can be accompanied by substantial ad- 
justments in the DM distribution, angular momentum, 
and other dynamic variables. The baryonic processes 
may also affect the survival of subhalos. This in turn 
can modify the basic parameters of the galactic disks 
which grow within DM halos, because disk-subhalo in- 
teractions can drive the disk evolution (e.g., Gauthier, 
Dubinski & Wilson 2006; Heller, Shlosman & Athanas- 
soula 2007a,b; Romano-Diaz et al. 2008b; Shlosman 2010 
and refs. therein). These effects are far from being fully 
quantified — a direct comparison between pure DM and 
DM+baryon models on subgalactic scales is limited both 
by numerical resolution and complex baryonic physics. 
In Paper I (Romano-Diaz et al. 2009), we compared the 
evolution of the prime halos of pure DM and DM+baryon 
numerical (hereafter PDM and BDM) models, evolved 
from identical initial conditions within the cosmological 
framework. Here we compare some aspects of the evolu- 
tion in the substructure associated with the prime halos 
in the PDM and BDM models. 

Recent efforts to understand galaxy formation have 
been spearheaded by pure DM halo formation in the 
cosmological framework, culminating in the large-scope 
Millenium (Springel et al. 2005), via Lactea (Diemand, 
Kuhlen & Madau 2007; Diemand et al. 2008), Aquar- 
ius (Springel et al. 2008) and Ghalo (Stadel et al. 2009) 
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simulations. Addition of a baryon component to the pure 
DM modeling is associated with difficulties related to a 
numerical modeling of dissipative processes, partly due 
to the sub-grid physics, and to our limited knowledge of 
physics of the ISM, star formation, energy and momen- 
tum feedback, etc., as shown in DM+hydrodynamic sim- 
ulations (e.g., Sommer-Larsen, Gdtz & Portinari 2003; 
Maccio et al. 2006; Berentzen & Shlosman 2006; Stin- 
son et al. 2006; Governato et al. 2007; Heller et al. 
2007a; Kaufmann et al. 2007; Romano-Diaz et al. 2009; 
Scannapieco et al. 2009), in chemodynamical simulations 
(e.g., Samland & Gerhard 2003; Brook et al. 2007a,b), 
or implementing semianalytical methods (e.g., Scanna- 
pieco et al. 2006). Comparison between various aspects 
of subhalo population between pure DM and baryonic 
models has been performed also for cluster scales (e.g., 
Nagai & Kravtsov 2005; Weinberg et al. 2008; Dolag et 
al. 2009). 

In Paper I and Romano-Diaz et al. (2008a, b) we quan- 
tified to what extent the baryons alter the DM density 
distribution within the prime halo. They contribute de- 
cisively to the evolution of its central region, leading ini- 
tially to an isothermal DM cusp, which is subsequently 
flattened to a DM density core — the result of heat- 
ing by dynamical friction of the DM+baryon subhalos 
during the quiescent evolution epoch. This confirmed 
previous work on this subject (e.g., El-Zant, Shlosman 
& Hoffman 2001; El-Zant et al. 2004; Tonini, Lapi & 
Salucci 2006; see also review by Primack 2009, as well as 
Johansson, Naab & Ostriker 2009). As a by-product of 
this process, the cold gas has been ablated from a grow- 
ing embedded disk, reducing the star formation rate by a 
factor of 10, and heating up the spheroidal gas and stel- 
lar components, triggering their expansion. We find that 
only a relatively small ~ 20% fraction of DM particles 
in PDM and BDM models are bound within the radius 
of maximal circular velocity in the halo, most of the DM 
particles perform larger radial excursions. We also find 
that the fraction of baryons within the halo virial ra- 
dius somewhat increases during the major mergers and 
decreases during the minor mergers. The net effect ap- 
pears to be negligible — an apparent result of our choice 
of feedback from stellar evolution. Furthermore, we find 
that the DM halos are only partially relaxed beyond their 
virialization. While the substructure is being tidally dis- 
rupted, mixing of its debris in the halo is not efficient and 
becomes even less so with redshift. The phase-space cor- 
relations (streamers) formed after z ~ 1 survive largely 
to the present time. 

This paper is structured as follows. Section 2 describes 
the initial conditions and the numerical modeling. Sec- 
tion 3 focuses on the basic properties of the PDM and 
BDM subhalos in our simulations, and section 4 com- 
pares the evolution of these subhalo populations in PDM 
and BDM prime halos. Discussion and conclusions are 
give in the last section. 

2. INITIAL CONDITIONS AND NUMERICAL MODELING 

For numerical details, model parameters and initial 
conditions, the reader is referred to Paper I. Here we 
only comment that the initial conditions for PDM and 
BDM models have been obtained using Constrained Re- 
alizations method (CRs, Bertschinger 1987; Hoffman & 
Ribak 1991; van de Weygaert & Bertschinger 1996). We 



followed the prescription of Hoffman & Ribak (1991) to 
build the initial conditions within a restricted cubic vol- 
ume of space with sides L = 8/i _1 Mpc in the ACDM 
cosmology, where a sphere of 5k" 1 Mpc is carved out 
and evolved from z = 120. Numerical simulations have 
been performed with the parallel version of FTM-4.5 
hybrid A-body/Smooth Particle Hydrodynamics (SPH) 
code (e.g., Heller & Shlosman 1994; Heller et al. 2007a; 
Paper I) using vacuum boundary conditions and physi- 
cal coordinates. The star formation (SF) modeling al- 
gorithm is described in Heller et al. (2007a) and in- 
volves multiple generations of stars, energy and momen- 
tum feedback from stellar evolution, supernovae and OB 
stellar winds. 

The total mass inside the computational sphere is 
~ 6.1 x 10 12 /i _1 M . To introduce the baryons, we have 
randomly replaced DM particles on the initial conditions 
grid by baryon (i.e., SPH) particles, so that Q m stays 
the same. The initial masses of DM and SPH particles 
are the same, 2.78 x 10 6 M . The evolution of various 
parameters characterizing the DM and baryons has been 
followed in 1000 snapshots, linearly spaced in the cosmo- 
logical expansion parameter a. 

2.1. Identifying and Measuring the Subhalos 

The HOP algorithm (Eisenstein & Hut 1998) used 
to identify the substructure within the computational 
sphere becomes inefficient deep within the prime halo, 
where density rises above the virial density p v ; r = 
A(z)p(z). Here A(z) is the critical overdensity at viri- 
alization (Bryan & Norman 1998) and p(z) is the back- 
ground density. The halo virial radius, i? v ir, is defined 
in the context of the spherical top-hat collapse model, 
M vir = 4/37rA(z) i o(z)i? vir 3 . Therefore, we divide the 
computational sphere into two regions — the boundary 
between these regions is i? v ir of a prime halo. In the 
outer, mostly unvirialized part we simply apply the HOP 
algorithm for halo finding. In the inner, mostly virialized 
part, we refine the method and perform the following it- 
erative procedure. 

• We start from R v „ and split the enclosed volume 
into arbitrary-shaped shells lying between isoden- 
sity surfaces, where p\+i = rjpi, with the shells 
counted from i? v j r inwards. We choose r\ = 1.5. 
Within each shell, the total (DM or DM+baryons) 
background density was taken as constant and 
equal to the average of two isodensities which sand- 
wich it. 

• The HOP algorithm was applied within each such 
shell, the subhalos have been identified and their 
tidal radii, R t , and masses, M t have been deter- 
mined from the condition that the subhalo density 
falls below the background (shell) density, p;. We 
check that such substructures are repeated in sev- 
eral iterations. Otherwise they are considered spu- 
rious and removed from the subhalo catalog. 

• In order to distinguish "real" structures from (un- 
bound) density enhancements, we require that real 
substructure will have a peak density p > r\p\ and 
that the number of DM particles associated with 
this candidate exceeds 100. We find that objects 
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Fig. 1. — Snapshots of PDM prime halo (top left) with the associated subhalo population (top right) and the DM in BDM (bottom left) 
prime halos with the associated subhalos (bottom right) at z = 2.33. The circle sizes scale with the tidal masses of the subhalos. The box 
size is 600 kpc. The colors correspond to the local DM mass density. M v ; r and R v i T of the main halos are ~2x 10 12 Mq and ~ 160 kpc, 
for PDM and BDM models. The tidal masses of three most massive subhalos are (8.8,4.3,3.7) X 10 10 M and (9.3,3.3,2.4) X 10 10 Mq, 
for PDM and BDM respectively. 



that include less than the above number of DM par- 
ticles provide unreliable statistical data in order to 
analyze their properties. 

• Finally, the background (i.e., unbound to a sub- 
halo) DM or baryon particles have been assessed 
using the velocity histogram of all particles within 
Rt in order to separate them from the genuine sub- 
halo mass contribution. We use velocities in the 
center of mass (CoM) of a subhalo. An alterna- 
tive method exists to separate the bound subhalo 
gas from the hot unbound background, using the 
bi-modal temperature distribution of the gas par- 
ticles. 

As stated above, we also define the virial and tidal radii 
of subhalos. For this we use the 3D density contours of 
subhalos obtained in the previously described procedure. 
The last contour still belonging to the subhalo is that 
whose density is defined in the first item. It is above the 
local background density of the prime halo (if located in- 
side the prime halo) or simply the background density (if 
outside the prime halo). The tidal radius of the subhalo, 
Rt, is defined as the maximal extent of such a contour 
from the subhalo CoM. The virial radius of a subhalo is 



defined similarly when the latter is situated outside the 
prime halo. The subhalo masses, both virial and tidal, 
are assumed to lie within density shells of these radii, re- 
spectively. Snapshots of prime halos and the associated 
subhalos are shown in Fig. 1 for PDM and BDM models 
at z = 2.33. 

3. RESULTS: BASIC PROPERTIES OF SUBHALOS IN PDM 
AND BDM MODELS 

We have compiled a comprehensive catalog of subha- 
los in PDM and BDM models which is used to infer their 
basic properties. In the following, we also choose repre- 
sentative samples of PDM and BDM subhalos out of the 
general catalog having different masses (i.e., in the mass 
ranges of 10 8 M© - 1O 9 M , 10 9 M Q - 1O 1O M , 10 10 M Q - 
10 11 Mq) at the time of crossing i? vu - of prime halos. Fi- 
nally, we look at the corresponding subhalo pairs in PDM 
and BDM. 

Our tests indicate that only subhalos with particle 
number exceeding 100 give a reliable estimate of their 
properties (see also Trenti et al. 2010). In order to iden- 
tify the time of a subhalo (tidal) dissolution, we mark 
the 20 densest DM particles in each of the subhalos and 
calculate the dispersion in their relative positions, cr(z), 
as a function of time (section 4). The tidal disruption of 
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Fig. 2.— Evolution of the SubHMF, N abh (M shh ), in PDM 
(blue, solid) and BDM (red, dashed) models using a power-law 



fit 



Power index a within R V1I is shown as 



N(M sbh ) - 

a function of redshift z and the cosmological expansion factor a. 
Only subhalos with N s bh > 100 have been used. This corresponds 
to M sbh £ 2.7 X 10 8 M e . The values of a at z = are -0.99 
and —1.13 for PDM and BDM respectively. We have smoothed the 
curves averaging them over the time frames. The averages of both 
curves are shown as a black dashed line at -0.86 (PDM) and a black 
dot-dashed one at -0.98 (BDM). 



a subhalo is reflected in a sharp increase in a. We choose 
a = 5 kpc for the time when the subhalo has been de- 
stroyed, and mark this time as z^es or Odes- This method 
has been tested and found as very reliable. 

Because we avoid fitting the NFW density profile to the 
subhalos, we modified the definition for the compactness 
parameter, c, used in the literature by replacing R s by 
Rvmax and i? v ir by i? t . This has been done by obtain- 
ing the rotation curves for the subhalos based on their 
total mass distribution. In Paper I we have shown that 
-Rvmax behaves largely similar to R s . The new defini- 
tion is c =Rt/Rvmax- Because i? V max is the most bound 
radius (i.e., maximum of the DM rotation curve), we typ- 
ically find that c <; 1 over the lifetime of a subhalo, i.e., 
it is destroyed before c drops below unity. 

3.1. Subhalo Mass Function 

We follow the evolution of the subhalo mass func- 
tion, iV(Af s bh) inside i? V i r of the prime halos, hereafter 
SubHMF. We approximate the SubHMF with a power 
law, N shh ~ Mg h (e.g., Moore et al. 1999; Ghigna et 
al. 2000; Springel et al. 2001),. We start by comparing 
the evolution of the power law index a in the SubHMF 
(Fig. 2). 

The main conclusions that emerge from Fig. 2 with 
regards to the evolution of a in PDM and BDM mod- 
els are as follows. First, the PDM SubHMF remains 
slightly shallower than the BDM one at all times, which 
is reflected in their averages as well, -0.86 vs -0.98 re- 
spectively. The final values at z = are -0.99 and -1.13. 
Second, the difference between these models, Aa(z), oth- 
erwise nearly constant, appears to increase and decrease 
after z ~ 0.5. During this time, the PDM SubHMF 
becomes slightly shallower, while the BDM one slightly 
steeper. Both functions converge thereafter. Third, the 
associated fluctuations in a{z) with respect to the mean 
values of a are related to the waves of inflowing low mass 
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Fig. 3.— Evolution of the SubHMF in PDM and BDM models: 
ratio of subhalo numbers (top) and masses (bottom) inside 30 kpc, 
within 30 kpc — R v in an d inside i? v i r of the prime halos. Only 
subhalos with JVdm > 100 have been used. All curves have been 
time-averaged over 50 frames. 



subhalos across i? v h- after the epoch of major mergers 
(Romano-Diaz et al. 2008a). Overall, this behavior is 
consistent with a ~ const, over the whole period, in the 
presence of an intrinsic scatter of ±15%. 

In principle, we cannot rule out that there is an evolu- 
tionary trend in a for z <> 0.2 in both models in Fig. 2. 
It may be also related to the overall depletion of the sub- 
halo population, i.e., the small number statistics, which 
would allow larger fluctuations from the mean toward 
the end of the simulations. As noted in the Discussion 
section, a similar trend was detected in the analysis of 
pure DM simulations of via Lactea (Madau et al. 2008). 

To compare the evolution of the subhalo populations 
within i? V j r and in the innermost 30 kpc in PDM and 
BDM, we display their number and mass ratios as a 
function of time in Fig. 3. Within the central 30 kpc, 
the BDM-to-PDM number ratio of subhalos stays above 
unity for the entire evolution after z ~ 3, while prior to 
z ~ 3 it oscillates in a narrow range around unity. After 
z ~ 1.5, this ratio is limited to 1.5 — 2.5 range and shows 
large amplitude fluctuations due to the clumpy influx 
(see above and in Romano-Diaz et al. 2008a), especially 
noticeable in the coresponding mass ratios in the lower 
frame of Fig. 3. On the other hand, this ratio remains 
steady and below unity in the less dense outer parts of 
the prime halos. This is also true for the subhalo popu- 
lation overall. 

The DM mass ratio, BDM/PDM, of subhalos (Fig. 3) 
within the inner 30 kpc appears to stabilize around unity 
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at later times, albeit fluctuating with substantial ampli- 
tudes. Beyond this region, i.e., between 30 kpc - -fiW, it 
stabilizes around 0.5 - 0.6, after the major merger epoch. 
This latter ratio also holds for the entire subhalo popula- 
tion, as we see from the black solid curve at i? v ir- Hence, 
there are more, by a factor of ~ 2, subhalos by number 
in the innermost region of the BDM prime halo, but not 
by mass. Their DM mass is comparable to that locked 
in the PDM subhalos. 

This excess of the subhalos within the inner 30 kpc of 
the prime BDM halo points to a shorter capture time 
of these objects in the central region than in the PDM 
model (as shown in section 4). This in turn leads to 
shorter lifetimes of the BDM subhalos. On the other 
hand, the total DM masses of PDM and BDM popu- 
lations in the central region are comparable. Partially 
this is because a fraction of DM, ~ 17%, has been con- 
verted into baryons in the initial conditions. As a result, 
a larger fraction of the BDM subhalos falls below the 
threshold of ./Vdm = 100, and, therefore, does not con- 
tribute to the subhalo mass in our estimates. On top 
of this, because the BDM subhalos concentrate deeper 
in the potential well of the prime halo, they are tidally 
stripped at a faster rate and hence have lesser masses on 
the average. 

3.2. Subhalo Population Spatial Distribution 

As a next step, we have computed the number den- 
sity of subhalos, n{R)/n v - l - Cl normalized by their aver- 
age density within i? v ir, in the PDM and BDM models 
(Fig. 4). As before, we only account for subhalos with a 
mass above the threshold mass of M s bh iS 2.7 x 10 8 M Q 
(i.e., -/Vdm > 100). The prominent feature in n{R) is a 
bump associated with some steepening of the slope at 
-Rvmax- As a next step, we analyze the DM mass distri- 
bution due to the subhalo population, p s bh, within the 
prime halos of both models. We find that the subhalo 
population distribution bears some similarity to that of 
the total DM mass distribution in the prime halos but 
also shows some differences. Within R vmaX} the DM den- 
sity decrease appears linear in the log p s bh — log R plot 
for both subhalo populations. We observe a sharp steep- 
ening of the density around R vmax . This characteristic 
radius is always larger in the PDM model by a factor 
of 2 (see Fig. 4 and Paper I). The subhalo population 
densities are typically steeper than -2 slope beyong this 
radius. We can trace the subhalos by almost a decade 
in radius inside i? V max, to the inner few kpc — the reso- 
lution limit of the subhalo detection. Within i? V max the 
slopes lie between -1 and -2, as seen in the snapshotes 
in Fig. 4 at various z, with the PDM slope being some- 
what steeper on the average (more about this below). As 
confirmed by the lower frame of Fig. 3, the rate of the 
subhalo influx into the central region is uneven which 
leads to an alternate dominance by the PDM or BDM 
subhalos at various times. Accordingly, the same behav- 
ior is found when we divide the subhalo population into 
spherical shells of equal width. 

We have fitted a time-averaged power law to the den- 
sity profile of subhalo population within R vmaX7 p s bh ~ 
r 7sbh shown in Fig. 5 (see Fig. 6). The PDM and BDM 
distributions differ, the former being somewhat steeper, 
with 7 s bh ~ — 1.5 vs — 1 for the BDM subhalos, confirm- 
ing Fig. 5. We note also that the profiles are more stable 



than those for the SubHMF, after the epoch of major 
mergers. The slopes are shallower than those of the to- 
tal DM distribution in the prime halos within i? V max and 
steeper outside these radii. This behavior can be domi- 
nated by dissolution of the subhalos in the central region. 

We have calculated the DM mass fraction of subhalo 
populations locked within R v [ r of the prime halos. After 
z ~ 3, both the PDM and BDM models behave similarly. 
This confirms our previous estimates of the DM associ- 
ated with subhalos and with their tidally disrupted re- 
manants (i.e., streamers, still visible in the phase space) 
within the prime halos (Fig. 11 of Paper I). In Paper I we 
have estimated that the total fraction of bound (to sub- 
halos) and unbound but not yet mixed material at z = 
to be ~ 8% — 9%. Here we find that about half of this 
DM is bound and the rest contributes to the streamers, 
i.e., to the not yet mixed material. 

We also comment on the axial ratios of the subhalo 
populations. As shown in Fig. 1, these ratios are much 
closer to unity in the presence of baryons. This agrees 
well with our results for the DM distribution in the prime 
halos in paper I. 

4. RESULTS: SUBHALO POPULATIONS IN PDM AND BDM 

MODELS 

Using the compiled catalog of subhalos in PDM and 
BDM models allows us to make a direct comparison be- 
tween corresponding pairs in order to understand the ef- 
fect of baryon presence on the subhalo population and 
its dynamics. Specifically, we have compared the corre- 
sponding subhalos trajectories, masses within their R t 
and -Rvmax, their concentration parameters, masses and 
mass ratios of all components to their total masses, spe- 
cific angular momenta, and determined their destruction 
times. The redshifts of the prime halo i? v i r crossing have 
been choosen to lie between z ~ 4 — 0.5. Two examples 
from this catalog are shown in Fig. 7 and provide a de- 
tailed quantitative analysis of these subhalo pairs. The 
subhalo population average properties are discussed later 
on. 

Analyzing the corresponding PDM and BDM pairs, 
we find that the redshift of a subhalo disruption, z^cs, 
is almost always higher for BDM than for the associ- 
ated PDM ones. This seems to be counter-intuitive be- 
cause one can argue that the BDM subhalos are more 
concentrated due to the baryon dissipation. This argu- 
ment, however, neglects the accompanying changes to 
the BDM prime halo, which becomes substantially more 
concentrated, as shown in Fig. 4 of Paper 1. The orbits 
of PDM and BDM subhalo pairs appear nearly identical 
up to their first pericenter. After this, the trajectories 
separate (Fig. 7a). The radial extent of the PDM trajec- 
tories remains larger, and the next apocenter of the PDM 
subhalo lies at larger R compared to the BDM. We also 
find that after the first pericenter is reached, the specific 
angular momentum of a PDM subhalo's particles is al- 
ways larger than that of the corresponding BDM, with 
respect to the CoM of the prime halo (Fig. 7d). In other 
words, the BDM subhalos settle on more radial orbits 
after the first pericenter passage. The tidal disruption of 
subhalos normally occurs after the pericenter, as can be 
seen from the suddent increase in o~(z) (Fig. 7b). The 
larger time spacing between the pericenters of the PDM 
subhalos is probably the main reason why the PDM pop- 
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Fig. 4. — Evolution of the subhalo population DM radial number density distributions within the prime halos in the PDM (blue solid) 
and BDM (red dashed) models, for subhalos with TV > 100 or M s t,h ^ 2.7 X 10 s Mq. The number density, n(R), is normalized by the 
average value of subhalo population DM density within R vlI . The prime halos DM density profiles are shown for comparison (same colors, 
dashed lines). Snapshots are shown for four different z. The vertical arrows correspond to the instantaneous values of fi Tmax (PDM and 
BDM) in these models with corresponding colors. All curves have been terminated at the instantaneous values of ii v j r for the prime halos. 



illation is characterized by a longer lifetime compared to 
their BDM counterparts. We return to this issue later 
on when we discuss the average properties of subhalo 
populations. 

The next important difference between the correspond- 
ing PDM and BDM subhalos is the evolution of the con- 
centration parameter. Already at higher z, close to the 
turnover radii, when the subhalos travel outside i? v ; r of 
the prime halo, c grows faster in BDM subhalos. This 
growth seems spectacular when compared to a nearly 
constant c of the PDM counterparts (Fig. 7c). It termi- 
nates typically when a subhalo enters i? V ir of the prime 
halo, or after the first pericenter of its orbit. However, 
the cutoff in the growth of c can come as early as the 
turnover radius of the prime. Frequently, the "normal" 
evolution of c is modified by interactions with other sub- 
halos. 

Before we turn to the average properties of the sub- 
halo populations, we note that the reason for consistently 
higher c in BDM subhalos is that their i? V max decreases 
early and substantially due to the accretion of baryons 
which concentrate deep within this radius (Fig. 7e), as 
can be confirmed by comparing the radial profiles of their 
circular velocities and by the maximal values of circular 
velocities (Fig. 7f). Hence, baryons clearly impact the 
internal structure of the BDM subhalos. We note that 



Rvmax in the PDM subhalos declines late in the evolution 
as well, typically in the last stage leading to a tidal dis- 
ruption (Fig. 7e) . This effect was first noticed by Klypin 
et al. (1999). However, c stays nearly constant in these 
objects, as the decline in i? V max is concurrent with the 
decline in R t . 

In order to follow the overall properties of subhalo pop- 
ulations in both models, we separate their evolution into 
4 spatial zones within and one zone outside the prime 
halos, namely, 0-0.25i? vir , 0.25 - 0.5i? vir , 0.5-0.75i? vir , 
0.75 — 1.0i? v i r and >i? v ir, where the virial radius is that of 
the prime halos (Fig. 8). In the innermost region, the av- 
erage c remains nearly constant in time for the PDM and 
slowly decreases for the BDM subhalos, below the PDM 
curve. For the next three regions, the average c is some- 
what more noisy for the PDM and much more noisy for 
the BDM ones. The overall trend is of a nearly constant c 
for PDM and a slight increase and a subsequent decrease 
for the BDM average values. However, the BDM curve 
always stays above the PDM one in the next three zones. 
This trend continues outside i? v ir, although the noise has 
disappeared completely and the separation between the 
average curves is larger. Strikingly, the dispersion of c 
values around the average is much larger for the BDM 
model in all zones, especially in the innermost one. 

As has been noted already in Fig. 7 for specific exam- 
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Fig. 5. — Evolution of the subhalo population DM radial mass distributions within the prime halos in the PDM (blue dashed) and BDM 
(red dashed) models, for subhalos with N > 100 or <; 2.7 X 10 s Mq. Similar but normalized number density distributions are shown 

in Fig. 4. Snapshots are shown for four different z. The prime halos DM mass distribution is given by the solid blue (PDM) and red (BDM) 
lines. Black solid lines show the -1 slope (upper) and -2 slope (lower) for a comparison. The vertical arrows correspond to the instantaneous 
values of i? V max (PDM and BDM) in these models with corresponding colors. All curves have been terminated at the instantaneous values 
of Rvir for the prime halos. 
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Fig. 6. — Evolution of subhalo population DM density profiles 
power index 7 s bhi within iJvmax shown in Fig. 5. The density 
profile is approximated by a power law, p s bh ~ r 7sbh for PDM 
(solid blue) and BDM (dashed red) models. The time averaging 
involves 20 frames. 



pies and confirmed in Fig. 8 for the whole subhalo popu- 
lation, the divergent behavior of the average c comes from 



a decrease in the average values of i?t and especially of 
Rvmax in the BDM subhalos. The maxima of the circular 
velocity curve move to much smaller radii in the BDM 
subhalos due to the captured baryons contribution. This 
is most dramatically demonstrated for the substructure 
outside the prime halo — here i? V max drops to a few kpc 
and stays nearly constant after the major merger epoch, 
while that of the PDM population is barely affected. 

To further compare the evolution of c for the PDM and 
BDM subhalos, we plot its values for individual subhalos 
at the time they enter i? v ; r of the prime halo and at the 
time they are tidally disrupted, c- ln and c<j es respectively, 
vs time (Fig. 9, left upper). The prevailing trends are 
obvious for both populations: no evolution for the PDM 
and a gradual increase by a factor of ~ 2 — 3 in BDM. 
Note that occasionally, the entering c exceeds that of the 
disruption, but overall they evolve in tandem. 

Because little evolution in c is observed in PDM, the c 
values lie in the corner of the available parameter space, 
Cdos — Qn, in PDM, but are much less confined in the 
BDM case (Fig. 9, upper right). Again we confirm that 
the average ratio of Cd es /cm is nearly constant for both 
populations. 

Comparison of the subhalo lifetimes following their en- 
trance into the prime halos shows that the average time 
before the tidal disruption increases in time during the 
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Fig. 7. — Examples of evolution of corresponding subhalo pairs in PDM and BDM models. Shown are (a) the radial trajectories given 
as the position of a subhalo with respect to the CoM of the prime halo, (fe) relative distances of the 20 densest particles, (c) concentration 
parameters, (d) specidic angular momenta with respect to the prime halo CoM, (e) i? vm ax, and (/) circular velocities at i? V max- Left frames: 
subhalo pair #1, Right frames: subhalo pair #11. 




Fig. 8. — Evolution of the concentration parameter c, Rt and Rvmax in PDM (blue) and BDM (red) subhalos within five radial zones: 
r = — 0.25-Rvir, T = 0.25 — 0.5ii v i r , r = 0.5 — 0.75-R v in r = 0.75 — li? v ir and T >Rvir (from left to right), where R vlI is that of the prime 
halo. The mean values for each population are shown as a yellow line (for population of red dots) and a light blue line (for population of 
blue dots). 



major merger epoch, then nearly levels off. The aver- 
age lifetime thereafter is about 1.5 Gyr for the PDM and 
~ 0.9 Gyr for the BDM subhalo population. However, 
the dispersion around this mean is substantial, especially 
during the major mergers epoch. For PDM, most of the 
subhalos are destroyed within 4 Gyrs, some survive much 
longer, ^5 — 9 Gyr (Fig. 9, lower left)). The lifetime of 



the BDM subhalos is much more limited, and most are 
destroyed within ^1 — 2 Gyr. 

Finally, we find that, on the average, the PDM sub- 
halos are destroyed within a narrow range in Cd cs , while 
the BDM ones exhibit a much larger dispersion in these 
values (Fig. 9, lower right). We also find a clear symme- 
try across the c m — Cdcs diagonal in both models — thus 




Fig. 9. — Evolution of the subhalos properties in PDM (left frames) and BDM (right frames). Upper left: c as a function of a; the 
average values of c at the time of entering the prime halo (blue solid) and being tidally disrupted (red solid) are shown as well. Upper 
right: concentration parameter at the time of the tidal disruption c^ cs vs concentration parameters at time of entrance, c ln ; the diagonal 
dashed line is shown for comparison. Lower left: lifetime of subhalos as a function of z and a; the average value are shown a solid blue 
line. Lower right: lifetime of subhalos as a function of c^ cs . Only subhalos that are destroyed before a = 1 are shown in these diagrams. 



the subhalos do not increase their c as they are tidally 
disrupted. This is interesting because Fig. 7c shows that 
c in BDM increase with time after penetrating the prime 
halo, although they fall sharply during the tidal disrup- 
tion process. 

Comparison with the PDM model reveals an important 
difference — although the BDM subhalos are more cen- 
trally concentrated than their PDM counterparts, they 
are disrupted at a higher rates. This applies of course 
only to the substructure that penetrates deep inside the 
prime halo, say the central ~ 30 kpc. 

5. DISCUSSION AND CONCLUSIONS 

The purpose of this work is to compare the properties 
and evolution of subhalos without baryons (PDM mod- 
els) and in their presence (BDM models). PDM and 
BDM models have been evolved from identical initial 
conditions designed by means of Constrained Realiza- 
tions method. The prime halos of both models have been 
compared in paper I. We have compiled a comprehensive 
catalog of PDM and BDM subhalos in order to analyze 
the statistical properties of these populations. In addi- 
tion, we created a subsample of corresponding PDM and 
BDM subhalo pairs to focus on the individual evolution. 
Lastly, we have performed some basic analysis of subhalo 
population orbits in PDM and BDM. 

Our main results are as following. We compile 
and compare the subhalo population mass functions 
(SubHMFs) in the range of z ~ 3 - 0. The SubHMFs are 
consistent with a single power law in the mass range of 
- 2 x 10 8 M - 2 x 10 11 Mq for both models. However, 



we note a systematic offset of the BDM population with 
respect to the PDM one. Namely, < a >~ -0.86 (PDM) 
and -0.98 (BDM), respectively. Some variablity in a(z) 
can come from an intrinsic scatter of ±15%. In the sim- 
ulations it is associated with waves of subhalos falling 
into the prime halos along the large-scale filaments. Such 
waves appear to be much more pronounced in the pres- 
ence of baryons. 

Next, we have analyzed the spatial DM number and 
mass densities of subhalo populations. We find that the 
latter one can be roughly described as a power law with 
a steepening around i? V max — the radius of maximal cir- 
cular velocity and hence largest binding energy. Inside 
-Rvmax the slope, 7sbh, is much more shallow than out- 
side this radius. A smoothed out 7 s bh is about -1.5 in 
the PDM and -1 in the BDM models, within i? vmax . The 
slope is closer to -3 outside this region. The shallower 
slope within i? V max is related to the elevated rates of sub- 
halo ablation and tidal disruption processes. Overall, the 
subhalo number density profiles follow the distributions 
of their respective main halos. 

Third, we find that the BDM population of subhalos is 
depleted at a higher rate compared to the PDM popula- 
tion. This result appears to be counter-intuitive because 
one expects the accreted baryons to increase the mass 
concentration of BDM subhalos and make them more 
resilient to tidal disruption. However, we find that the 
trend leading to the increased mass concentration of the 
prime halo in the BDM model dominates and ultimately 
contributes to the dissolution of subhalos. Specifically, 
trajectories of PDM-BDM subhalo pairs appear to di- 
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verge after the first pericenter. The subsequent BDM 
subhalos trajectories become confined to the innermost 
regions of the prime halo where they experience increased 
dynamical friction and spiral in towards their ultimate 
tidal disruption. The lower slope for the inner density 
profile of the BDM subhalos compared to the PDM stems 
from their elevated rate of the tidal disruption. The small 
number of central subhalos leads to significant fluctua- 
tions in the slope. Yet, on the average, the BDM popula- 
tion is depleted faster in this region as shown by Figures 
6 and 9, resulting in a lower slope. 

We now compare our results with those found in the 
literature, starting with the SubHMFs. Most of the mod- 
eling of SubHMFs on the galactic scales has been focused 
on the pure DM models. Comparison with the BDM 
models is nearly always limited to galaxy clusters. In the 
galaxy and cluster mass range, the slope a was found to 
be consistent with -1, for the PDM models (Moore et al. 
1999; Ghigna et al. 2000; Springel et al 2001; Stoehr et 
al. 2003; De Lucia et al. 2004; Gao et al. 2004; van den 
Bosch et al. 2005; and Diemand et al. 2007, 2008). More 
recent multi-resolution Aquarius simulation (Springel et 
al. 2008) estimate the range for a to lie between -0.87 
and -0.93, close to our value of -0.86. 

The later evolution of the PDM SubHMFs has also 
been noticed in the Via Lactea simulation. For via Lactea 
and additional simulations, Madau et al. (2008) showed 
that the slope of the mass function at two different red- 
shifts (z — 0.5, — the only ones published) has slightly 
changed, steepening from —0.92 to —0.97. Such a "last 
moment" steepening is compatible with our Fig. 2, as 
noted in section 3.1. We find that this possible steep- 
ening in PDM and BDM subhalo populations can be re- 
lated to the small number statistics toward the end of 
the simulation — i.e., we possibly observe a large am- 
plitude variation around the mean for the slope a. We, 
therefore, leave this issue open. 

There is much less data on the SubHMF in the presence 
of baryons, and it deals with the mass range of clusters 
of galaxies. Comparing simulations (with and without 
baryions) drawn from the same initial conditions, Wein- 
berg et al. (2008) found the two subhalo mass func- 
tions to be very similar. They conclude that the dissi- 
pative baryonic component has only a "small" impact 
on this global measure of the subhalo population. On 
the other hand, Dolag et al. (2009) find that the ad- 
dition of baryons somewhat modifies the SubHMF for 
galaxy clusters, with cspdm to be somewhat steeper than 
c^bdm- While their values of a are similar to ours, the 
order appears inverted. Finally, for z = 0, Libeskind et 
al. (2010) point out that the BDM substructure is more 
radially concentrated than in the PDM. 

Next, we turn to the DM density profiles of the sub- 
halo population. We find Maccio et al. (2006) as the 
only work on this issue within the galaxy mass range, 
although they stopped the gas cooling at z — 1.5 (see 
also Sales et al. 2007, although they only analyze the 
most massive subhalos within their simulations). This 
work provides only the subhalo number density profile 
(their Fig. 6), which is consistent with ours. Similar re- 
sults have been obtained by Nagai & Kravtsov (2005) and 
by Weinberg et al. (2008) — all deal with the subhalo 
population number density profiles, albeit in the galaxy 
clusters mass range. Furthermore, we find that there is 



little spatial bias (less than 10%) between the subhalo 
population DM mass density and the DM distribution of 
the main halos (e.g., Sales et al. 2007). 

Extending the subhalo population number density to 
the DM mass density, we find that a new feature — steep- 
ening of the mass density profile, appears to be associated 
with the position of R vmax in each model. Around i? V max, 
the slopes change and steepen from 7 s bh ~ [— 1] — [— 1-5] 
to below -2. If verified in other models, this is an interest- 
ing effect related to increased ablation and dissolution of 
subhalos in the cnetral region of the prime halos. It has 
implications on dynamics of subhalos and reflects their 
trapping within i? V max- 

Turning to the properties of individual subhalos in 
PDM and BDM models, specifically to the evolution of 
the concentration parameter, c, we find that except in 
the innermost 30 kpc the average c of the BDM popula- 
tion is higher than that for the PDM one (Fig. 8). Hence, 
it seems surprising that both PDM and BDM subhalos 
have the same c when they enter the prime halo and when 
they are tidally disrupted, i.e., Cd cs ~ Q n , as shown by 
Fig. 9. Resolution of this puzzle is provided by Fig. 7c. 
Indeed, c in BDM stays well above that for PDM and 
increases with time. However, in the process of a tidal 
disruption, it drops sharply, and in many cases to near 
the original value, thus assuring that the population av- 
erage has Cdcs ~ Qn- This coincidence is related to our 
choice of the "moment of disruption" - taking a — 5 kpc 
for the dispersion in the positions of 20 densest particles 
in a subhalo. However, this definition has been tested, 
including visual test as well as a test in the phase space, 
and found highly reliable. Any other choice of the critical 
value for a would be more ad hoc. 

Our result for the average c for the subhalos population 
in both models agree well with the conclusion by Dolag 
et al. (2009) that BDM population is on the average 
more concentrated than the PDM, albeit they show only 
a single snapshot in redshift and in the cluster regime. 
They also use the standard definition of c by fitting the 
NFW profile to subhalos. However, we find that the 
subhalos are tidally truncated and avoid using the NFW 
model fit. The tidal effects are especially visible in the 
central region of 30 kpc where we find that the PDM arc 
on the average less concentrated than the BDM ones. 
This of course must be taken together with the caveat 
that the PDM subhalos spend less time inside this region, 
as discussed above (see also Fig. 7a). 

This brings us to an interesting point, namely, which 
subhalo population is tidally disrupted faster. While it is 
obvious that baryons provide an additional "glue" to the 
DM subhalos, they also affect the innermost region of the 
prime halo. This is obvious from our analysis of DM dis- 
tribution in the prime halo of the BDM model (Paper I), 
where i? V max is essentially halved compared to the PDM 
prime. Even more prominent changes are affecting the 
cusp region (Romano-Diaz et al. 2008a). Clearly, the res- 
olution of the question of which population is destroyed 
at a higher rate depends on the competition between two 
effects: increase in the concentration of the subhalos ver- 
sus increase in the concentration of the prime halos. In 
our simulations, the BDM population appears to be de- 
stroyed at a higher rate in the central region. However, 
if the energy and momentum feedback from stellar pop- 
ulation and from the central supermassive black hole are 
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able to decrease and maintain a lower fraction of baryons 
in the prime halo, the effect on the DM density profile will 
be sharply decreased and the BDM subhalo population 
will be much more resilient to destruction than the PDM 
one. It is possible that nature benefits from both solu- 
tions. In addition, it is plausible that the galaxy cluster 
environment prefers the survival of the BDM subhalos, 
(as in Weinberg et al. 2008; Dolag et al 2009), while on 
the smaller mass scale, the galaxy environment lead to 
the opposite. 

With respect to the total fraction of the halo mass 
invested within substructure we find that ~ 5% is locked 
within subhalos. This result is consistent with previous 
published analysis (e.g., Ghigna et al. 2000; Springel 
et al. 2001; Stoehr et al. 2003) which provide a range 
between 5% and 20%. However, there is no agreement 
on this matter, as Moore et al. (2001) argues that the 
true fraction might approach unity if the subhalos could 
be identified down to very small masses. Our results are 
also consistent with the Via Lactea simulations where the 
subhalo mass fraction is ~ 5.3% within i? v i r , but appear 
somewhat lower compared to Aquarius simulation, which 
quotes ~ 11.2%. 

To summarize, we have analyzed the properties and 
evolution of subhalo population in the mass range of 
10 s M Q - 10 11 M Q in pure DM and DM+baryons mod- 
els evolved from identical initial conditions. Our main 
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results are that the subhalo mass functions can be fitted 
with a single power law which is compatible with no red- 
shift evolution between z ~ 3 — 0, but which is somewhat 
steeper in the presence of baryons. We have also com- 
puted the number and mass density profiles for the DM 
subhalo component in both models. We find that the 
DM mass density is shallower inside the radius of the 
maximal circular velocity in the prime halos and steep- 
ens outside. Finally, we compared the disruption rates 
for the subhalo populations and find that, in the pres- 
ence of baryons, the subhalos are destroyed at a higher 
rate within the central 30 kpc of the prime halo — a di- 
rect result from the increased central mass concentration. 
However, with a larger feedback from stellar populations 
and the central supermassive black holes, the effect can 
be reversed. 
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